NASA Technical Memorandum 102344 
TCOMP-89-22 


A New Uniformly Valid Asymptotic 
' Integration Algorithm for Elasto-Plastic- 
Crcep and Unified Viscoplastic Theories 
Including Continuum Damage 


_S 


A. Chulya 

Institute for Computational Mechanics in Propulsion 


: Lewis_ Research Center 
Cleveland, Ohio 




and __ _ 

K.P. Walker 

Engineering Science Software, Inc. 
Smithfield, Rhode Island 


December 1989 






E - 5052 


A NEW UNIFORMLY VALID ASYMPTOTIC INTEGRATION ALGORITHM 
FOR ELASTO-PLASTIC-CREEP AND UNIFIED VISCOPLASTIC 
THEORIES INCLUDING CONTINUUM DAMAGE 
A. Chulya 

Institute for Computational Mechanics in Propulsion 
NASA Lewis Research Center 
Cleveland, Ohio 

and 

K.P. Walker 

Engineering Science Software, Inc. 

Smithfield, Rhode Island 

ABSTRACT 

A new scheme to integrate a system of stiff differential equations for 
both the el asto-plasti c-creep and the unified viscoplastic theories is 
presented. The method has high stability, allows large time increments, and 
is implicit and iterative. It is suitable for use with continuum damage 
theories. The scheme was incorporated into MARC, a commercial finite element 
code through a user subroutine called HYPELA. Results from numerical problems 
under complex loading histories are presented for both small and large scale 
analysis. To demonstrate the scheme's accuracy and efficiency, comparisons to 
a self-adaptive forward Euler method are made. 

NOMENCLATURE 

[B] strain-displacement transformation matrix 

[C3 stress-strain material property matrix 

E Young's modulus 

e.j deviatoric strain tensor 

ePj deviatoric plastic strain tensor 

e c effective creep strain 

e p effective plastic strain 
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Supers 


drag stress 

global stiffness matrix 
Inelastic strain increment 

i. u 

vector of unbalanced force at i tn iteration 
apparent area 
deviatoric stress tensor 
net area 

vector of increments in nodal point displacements at 
vol ume 

vector of stress 
Kronecker delta 
total strain tensor 

creep strain tensor 

elastic strain tensor 

plastic strain tensor 

inelastic strain Incremental vector 

Poisson's ratio 

von Mises effective stress 

Cauchy stress tensor 

net stress 

Instantaneous yield stress 
damage parameter 
back stress tensor 
cripts: 


T 


transpose 



o 


Initial value 


• rate 

INTRODUCTION 

The Increasing demand for Integrity and reliability In metallic structural 
components that have been subjected to a complex cyclic thermomechanical 
environment has stimulated the improvement of Inelastic analysis methodology. 
The aerospace Industry very much needs durable and efficient combustor and 
turbine structural components In the modern gas turbine engine. Creep and 
fatigue cracking and creep buckling distortion of combustor liners caused by 
high temperatures and impact and erosion damage decrease turbine durability; 
this has led to the development of a phenomenological theory of unified 
constitutive equations that describe time and temperature dependence in the 
plastic regime, in contrast to the time-independence of classical plasticity. 

In the past the inelastic strain comprised a time-independent (plasticity) and 
a time-dependent (creep) term; these terms were calculated by using classical 
plasticity and creep theories, respectively. However the physical Interaction 
between creep and plasticity was observed through several deformation 
phenomena, that is, cyclic hardening or softening, creep recovery, and rate 
sensitivity. The unified constitutive theory Is considered to be superior in 
predicting and governing the physical process, as compared to the classical 
plasticity-creep theory. 

Although neither theory has been widely applied in structural analysis of 
samples under complex loading histories, the unified constitutive theory has 
been especially neglected. This is due mainly to difficulties associated with 
the system of very stiff differential equations in certain regimes. This 
mathematical stiffness requires use of a very small time step In order to 
integrate the constitutve models without loss of stability. As a result. 
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computation time becomes enormous, and under complex loading, solving the 
problems often becomes impossible. 

The Importance of an efficient aloglrthm to Integrate the Inelastic 
constitutive models Is obvious. An explicit algorithm, such as a forward 
Euler In conjunction with a self-adaptive scheme, has been widely used largely 
because it is simple and the computation is inexpensive. However, this 
algorithm Is a subi ncremental type, which is noniterative In nature. In this 
case, convergence of the solutions depends significantly on the judgment of 
engineers, who tend to be conservative. An implicit alogorithm, which is an 
iterative type, is more stable amd accurate but prohibitively expensive. In 
recent years, several approaches have been developed to make the algorithm 
less dependent on the analysts. Banthia and Mukherjee (1982), with their 
one-step Euler integration scheme with a variable time step. Improved the 
scheme by Imposing a better time-step control. This approach takes advantage 
of the fact that the equations appear to be stiffer for large strain rates. 
They chose a time step that gives more accurate results, but it Is slightly 
less efficient than their previous algorithm. Miller and Tanaka (1988) 
developed a noniterative, self-correcting solution (NONSS). Their method is 
similar to the Newmark 3-method in that a parameter that determines whether 
the method is explicit or implicit Is introduced. This method reduces to the 
forward Euler method when 3 > 0. Implicit quantities are removed in the 
NONSS method by Taylor expansions of state variables. The NONSS method is 
unconditionally stable of 32 1/2, but it requires setting up a Jacobian 
matrix and solving a set of linear equations at each time step. Accuracy is 
maintained through self-adaptive time control and by correcting errors at the 
current step. Since this method has been used in one-element applications 
only. Its applicability to finite element analysis remains to be seen. 
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Despite these efforts, a generally applicable method that is implicit, 
iterative, stable, and inexpensive as well as convenient for implementation 
into finite element codes has not yet been developed. The objective of this 
report is to demonstrate such a method. The proposed method is based on 
transforming the differential equations of constitutive models to an Integrated 
form as proposed by Walker (1976, 1980). These integrated equations are then 
approximated by uniformly valid asymptotic expansions (UVAE). A concise 
mathematical derivation is presented for both the classical theory of 
plasticity and creep and the unified viscoplastic theory. The advantage of 
this method in continuum damage mechanics is presented as well. Implementation 
into the commercial MARC finite element code is demonstrated. Results of 
numerical examples for small- and large-scale problems at high temperatures are 
then presented. Comparisons to the self-adaptive forward Euler (SAFE) scheme 
are made as well, to show the accuracy and efficiency of the proposed 
Integration scheme. 

DIFFERENTIAL FORMS OF ELASTO-PLASTIC-CREEP AND UNIFIED VISCOPLASTIC 
CONSTITUTIVE EQUATIONS 

A fundamental observation when comparing elastic and inelastic analysis 
is that for elastic solutions the total stress can be determined from the total 
strain alone, whereas in an Inelastic response calculation the total stress 
beyond the yield point depends on both the stress and strain histories. 

Typical inelastic phenomena are plasticity, creep, and viscoplasticity, and a 
very large number of material models have been developed in order to 
characterize such material response. In this section the basic forms of 
differential equations for both the classical theory of plasticity and creep 
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and the unified viscoplastic theory are presented. A brief review of the well- 
known, self-adaptive forward Euler (SAFE) Integration algorithm Is included In 
the last part of this section. 

Classical Theory of Creep and Plasticity 

The simplest and most widely used material model by far employs the 
classical plasticity theory to characterize short-term deformation and the 
classical creep theory to characterize long-term deformation. The differential 
equations are formulated explicitly and independently. For small displacement 
and small strain formulation, the total strain rate e.. Is decomposed into 

J 

elastic, plastic, and creep strain rates, 


•e 

ij " e ij 
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where = component of total strain rate tensor, = component of elastic 
strain rate tensor, = component of plastic strain rate tensor, 

= component of creep strain rate tensor. The constitutive law for an 
Isotropic material with temperature-dependent moduli (Fung, 1965; Malvern, 
1969) Is 


where 
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C. . = X S. . S + i/s, S. + S, 5. ^ (3) 

i J r s ij rs 1 r js Is jr J 

\ = Ev/(1 + v)(l - 2v); p = E/2(l + v) ; E = Young's modulus; v = Poisson's 

ratio; and 6.^ = Kronecker delta. 

The plastic strain rate is calculated by using the classical theory of 
time-independent plasticity (Meldelson, 1968; Fung, 1965; Malvern, 1969). The 
von Mlses yield function for noni sothermal , isotropic haradening can be written 
as 
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where Is the instantaneous yield stress and 
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e p is effective plastic strain, 
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elj’j is deviatoric plastic strain, and 0 is the temperature. With the yield 
stress defined as a function of the effective plastic strain eP and 
temperature 0, Eq. (6) can be used directly to evaluate f (Snyder, Bathe, 


1981). 

The creep strain rate Is determined by using a modified equation-of-state 
approach. This approach Includes strain hardening for variable loading and the 
Oak Ridge National Laboratory's auxiliary hardening rules of cyclic behavior 
(Pugh et al., 1972). The final result is written as 


e 1j * rS ij 

where T is a scalar variable, 



a is von Mises effective stress. 
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e c Is the effective creep strain rate, 


m 

e c = A Q (o) ] (t) 


(v) 


(id 


t is the effective time, which can be determined iteratively, and Aq, m] , 
and m 2 are given constants. Equation (11) is derived from the uniaxial 
creep law, which has been generalized to multiaxial conditions by utilizing 
the effective stress and effective creep strain. Other creep laws could be 
employed in a similar manner. 

The final inelastic strain rate can be written as 


-1 

: U 


• p *c 

e 1j + e 1j 


or 


( 12 ) 


■ <ft ns u ; 

Unified Viscoplastic Theory 

The new theories to characterize material behavior at high temperatures 
are known as unified theories, in the sense that plastic and creep strains are 
considered as arising from the same physical mechanism. One or more state 
variables are Introduced in the constitutive equations. Again a small 
displacement and small strain formulation is used, with the total strain rate 
decomposed into elastic and inelastic ejj parts. 


: 1j 


• e *1 

e ij + e 1j 


(13) 


The relation between the elastic strain rate and the stress rate Is given 

Hooke's law (see Eqs. (2) and (3)). 

The general form of the unified viscoplastic constitutive equations can be 
written as, 
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where Q.^ and K are the new state variables. The Is defined as the 

back (equilibrium) stress tensor, and K Is the drag stress. The R and B 

are written as 

1/2 


' R * G c 1j e ij) 

= [I ( S 1j ' ! Q 1j)( S ij " f Q ij)] 
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The f n Is a function of B/K to the power n, where A^ , A^, and n are 

• • 

material parameters, A^ is the strain hardening function, and G and J are 
recovery functions. 

Equation (14) is the flow law defining Inelastic strain rate as a function 
of applied stress, state variables, and temperature; this function was selected 
to represent creep curves. As a result, three extensively used functions in 
creep theories namely power, exponential, and hyperbolic sine functions, are 
adopted in the unified viscoplastic theories. Of these, the power function has 
been broadly used because of its numerical simplicity. 

Equations (15) and (16) are known as evolutionaly equations and are 
generally written in the context of a hardening and recovery form. The strain 
hardening function varies according to the value of the inelastic strain rate. 
The recovery function can be divided into dynamic and static recovery 
components . 
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In this report, the specific constitutive equations proposed by Walker 
(1976, 1980, 1981), Krl eng-Swearengen-Rohde (1978), and Miller (1976) were used 
to test the new Integration scheme; their equations are given in detail in 
Appendixes A, B, and C, respectively. 

Self-Adaptive Forward Euler Integration Algorithm 

Since the constitutive equations presented in the previous sections 
represent a system of first-order nonlinear differential equations, they can 
be written in a compact form, with the assumption that the stress and strain 
field are known at a given time t, as 

y = f(y,t) (19) 

where y represents the vectors of stress (Eq. (2)), inelastic strain 
(Eq. (12) or (14)), back stress (Eq. (15)), drag stress (Eq. (16)); and f(y,t) 
is an abbreviation for the nonlinear functions on the right side of these 
constitutive equations. These differential equations are, in general, highly 
nonlinear and have stiff regimes, particularly in the viscoplastic theory. 
Consequently, a very small time step is often required in order to use standard 
numerical Integration techniques to solve these equations without loss of 
stabi 1 i ty. 

Various numerical Integration methods have been proposed for solving 
siffness problems. Most of them are, however, Intended for fields other than 
structural mechanics. For example, Gear (1971) developed a program for 
handling a general class of stiffness. Although Gear's methods have been 
successfully appilied to uniaxial one-element analysis (Miller, 1975), his 
package is not suitable for a large-scale finite element analysis because of 
its extensive computer time and storage requirement. 

In the context of finite element analysis of rate-related problems, a 
number of numerical methods have been recommended. Numerical comparisons were 



extensively performed by Chang (1985) and Lindholm et al. (1985a, 1985b). The 
SAFE method has been found to be computationally efficient, especially when 
connected with the sub 1 ncremental approach (Cassentl , 1983a, 1983b). The most 
significant advantage of this approach Is that the numerical Instability 
Incurred in using an explicit method has been diminished in such a way that no 
reduction in gobal step size Is necessary. The local step size is adjusted on 
the basis of a comparison between an estimated error and prescribed error 
bounds. This scheme Is used for comparison purposes In the section Numerical 
Examples and Comparisons. The details of the scheme follow. 

For the solution of Eq. (19), the Initial values y(t = 0) = y^ have to 
be prescribed. The numerical solution is performed in discrete time steps 
At, starting from a known solution at time t. This time step At Is the 
current finite element global load Increment and Is divided into NSPLIT equal 
subincrements. The Integration of the system In Eq. (19) is then performed by 
using forward differences with a smaller step size as 

y t+At = y t + (NSPLn) f<t> (20 ‘ 


The above equation is repeated NSPLIT times and the solution of y at time 
t + At is obtained. There are three possible ways to determine NSPLIT, 
depending on the magnitude of the change In a strain measure for every 
subincrement. The change In the strain measure Is defined as 


where 
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If the value of ERROR is between the specified tolerances, usually 1 x 1 0~ 4 and 
lxlO“ 5 (defined as ERROR 1 In the section on numerical examples), then there Is 
no change In NSPLIT for the next subincrement. However If ERROR is less than 
the upper bound tolerance, NSPLIT Is reduced by half. When ERROR is greater 
than the lower bound tolerance, NSPLIT Is doubled and the current subincrement 
step Is repeated. If NSPLIT exceeds the maximum number specified, the 
conventional explicit forward Euler scheme Is exploited with a fixed number of 
subincrements throughout. This scheme, often called "successive substitution" 
requires that a very small time step be enforced in stiff regions to avoid 
numerical instability. Note that this Is not an Iterative scheme. 

INTEGRAL FORMS OF ELASTO-PLASTIC-CREEP AND UNIFIED VISCOPLASTIC CONSTITUTIVE 
EQUATIONS 

The fundamental concept in deriving a uniformly valid asymptotic 
integration scheme is to convert the constitutive differential equations 
presented in the preceding section into Integral form. In this section the 
procedure for transforming a differential form into an integral form is 
presented for el asto-pl asti c-creep and unified viscoplastic theories. 
Elasto-Plastl c-Creep Constitutive Integrated Equations 

The inelastic strain rate tensor is written in the form of the deviatoric 
strain rate tensor e.^ as 
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Equating Eqs. (24) and (27) yields 
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(27) 


(28) 


2pe 1j " S ij = QS 1j 

Rearrange the above equation to give a form of a first-order differential 
equation 


(29) 


s tj + QS tj ■ 2| *ij 
Integrating Eq. (30) for S^ at time t yields 

S.j(t) = S. j ( 0 )expt— Q< t ) ] + J t exp - J* dx 


(30) 


3 e i . 

2 M - gl 1 d ^ 


(31) 


Since at t = 0, Sjj(O) = 0; and e^ = - S }j e kk/3; then Eq. (31) becomes 

the final Integral form of Eq. (27) 
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Equation (32) represents the Integral form of the differential equation 
defined by Eq. (12), and it has a new scalar parameter Q. 

Unified Viscoplastic Constitutive Integrated Equations 

Each state variable In the differential form of viscoplastic relations 
presented In the section Unified Vfscoplastlc Theory Is converted into the 
integral form by using a procedure similar to that for elasto-plastic creep. 
Miller's model with three state variables demonstrates these transformations. 
Inelastic strain. - The Inelastic strain rate tensor is 
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where 

B Isa material constant, and 0' Is defined In Appendix C. Equation (35) 
can be rewritten as 


where 
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The second Invariant of the Inelastic strain rate tensor Is written as 


(37) 
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Substitution of Eq. 


S-G-Wi) 

(37) into Eq. (39) gives 
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From Eq. (40) the relation Z = 3yR/L may be substituted into the denominator 
of Eq. (38) to give 


from which the relation 



2h 2 m 



R = 
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is obtained. 

Notice that Eq. (41) can be extracted directly from the right side of 
Eq. (35). Rearranging Eq. (41) gives 


Z = 


K 
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(42) 


which, when substituted into Eq. (38), yields 
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Equating Eq. (37) with Eq. (24) forms 
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y 1j = S ij " 3 
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Then Eq. (44) becomes 


y fj * Ly,j - 2 v i^ - § Q, j 


(46) 


Equation (46) is in the form of a first-order differential equation and can be 
integrated as 

,t I pt 
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where yij(O) is the initial value. If Sij(O) = 0 and Qjj(0> = 0, then 
yij(0) = 0. By substituting Eq. (45) and the properties of deviatoric stress 
and strain into Eq. (47), it can be written as 
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Back stress. - The back stress in differential form Is 
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The symbols n, Hp A 1 , B, Q* , and k are material constants independent of 
temperature; T is the temperature In Kelvin, and T^ is the melting 
temperature of the material. We can rewrite Eq. (50) Into a first-order 
differential equation as 

•1 
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Equation (51) can be integrated to give 
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where Q^(t) » 0 at t = 0, and 
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Drag stress . - The final state variable of Miller's theory is the drag 
stress and is defined In differential form as 
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K = H 2 R 


1/2 A 

C 2*(§Vu) ■ 57 1(3 - H 2 C 2 Be'[ S lnh(A 2 K 3 )] 


(55) 


where H 2 , C 2> and A., are material constants and independent of temperature. 
Again Eq. (55) can be rearranged into the first-order differential equation as 
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Integrating Eq. (56) under the Initial condition K(t) = Kq at t = 0 gives 
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A UNIFORMLY VALID ASYMPTOTIC EXPANSION INTEGRATION ALGORITHM 

With the Integrated form of each state variable presented in the previous 
section, the asymptotic expansion can now be used to represent each integral 
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of both the elastic-plastic creep and the unified viscoplastic theories. The 
final forms appear in recursive relations and need to be solved by a Newton- 
Raphson iterative technique. 

El asti c-Pl asti c— Creep Uniformly Valid Asymptotic Expansion 
Equation (32) written at time t + At can be shown as 

a..j(t + At) = (x + | p)s i je kk <t + At) + Ijj<t + At) (60) 
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Because of incremental formulation in nonlinear analysis, Eq. (61) can be 
separated into two parts as 


I i j ( t + At) = I!j(t) + I j j (At ) 

where 
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The unity expression ^e^^e is Introduced into Eq. (63) to give 
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From Eqs. (60) and (61), the right side of Eq. (65) can be identified as 
I M (t); thus 

Ijj(t) = exp(-AQ)I. ^(t) 


or 




= exp(-AQ) 


L'U (t> 


- ( x * f ^) s i j c kk ct> ] 


where 


( 66 ) 


AQ = Q(t + At) - Q(t) = Q(t + At)At (67) 

The only integral in Eq. (62) is now I!^(At). According to Eq. (A6) in 
Appendix A (Walker, 1987), this integral can be represented by a uniformly 
valid asymptotic form. If only the first term of Eq. (A6) is used, the 
approximated recursive relation of Eq. (64) is 

- (2m ^ - f M as, Kk )[ ’ - e ^ < - aQ> ] (68) 
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The asymptotic recursive form of Eq. (60) becomes 
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Unified Viscoplastic Uniformly Valid Asymptotic Expansion 


A procedure similar to that for the elasto-plastlc-creep model Is used to 
obtain the final UVAE recursive forms of viscoplastic model for all three state 
variables, namely, Cauchy stress, back stress, and drag stress. 

For Cauchy stress, the relationship Is 

a. j ( t + At) = | Q.j(t + At) + (x + | p)s i; .e kk (t + At) 
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For back stress, the relationship Is 
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(74) 
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Finally for drag stress, the relationship is 
K(t + At) = < 0 + [k( t) - K q ] exp(-AJ) + H 2 <jc 2 + [| Q^(t + At)Q, j<t + At) 


, 1/2 


p K 3 < t^ AR [— 


- e ;p ( - - A ^l (75) 
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where 


AJ 


At 


K(t) - K Q 


H 2 C 2 Be '{ sinh [ A 2 K3(t> ] / 


(76) 


AR = R(t + At)At 


Newson-Raphson Iteration 

Unlike the forward Euler Integration scheme, Eqs. (70), (71), (73), and 
(75) are recursive in nature. Each unknown state variable at time t + At 
Involves a single parameter (i.e., AQ, AG, or AJ) which. In turn, requires a 
knowledge of the parameter's unknown state variable. These equations are the 
recursive or implicit equations. Therefore a technique such as the Newton- 
Raphson iteration is required. However this Newton-Raphson implicit iterative 
scheme is different from the implicit integration scheme of differential 
equations (Chang, 1985) in that its Jacobian matrix is much smaller. For the 
e 1 as to-p 1 as 1 1 c-creep theory, instead of iterating over six components of Cauchy 
stress, Eq. (70) iterates over one parameter, AQ. In the case of unified 
viscoplastic theory, the size of the Jacobian matrix is reduced from 13 x 13 
(six components of Cauchy stress, six of back stress, and one of drag stress) 
to 3 x 3 (AQ, AG, and AJ); this is the tremendous advantage of transforming 
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differential equations to UVAE equations. Since the Jacobian matrix of the 
elasto-plastlc-creep model is a subset of that of the viscoplastic model, the 
latter will be used to demonstrate the Newton-Raphson technique. 

The governing Iterative equations are obtained from Eqs. (72), (74), and 
(76) for each state variable and written in the function forms as 


f,<aO, ag. aj> , AO _ 


At 


1 / D \ 1 

slnh" < — 


VB0 
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2/3 


f 2 ( AQ , AG, AJ) = AG - H ] B9‘ ^si nh ^ £)„ 
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n n 
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(t + At)G..(t + At) 

J J * 




(77) 


f,(AQ, AG, AJ) = AJ - t — T [ H-C-B0 

j I IX / \ IS \ \ £ 
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The iteration starts with judiciously chosen initial guesses for AQ, AG, and 
AJ. The intent is to find the solution of the equations 


f m (AQ, AG, AJ) = 0 


or (78) 

f mK) * 0 

where m = 1,2,3 and AU^ is a vector of AQ, AG, and AJ. Equation (78) 
represents a system of nonlinear equations. The most frequently used iteration 
scheme for the solution of these equations is some form of a Newton-Raphson 
iteration (Stricklin et al . , 1973; Oden, 1972; Bergan et al . , 1978). By using 
a Taylor series expansion and retainning only the first-order term, the 
iterative form of Eq. (78) is written as 
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(79) 
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(80) 


and 



(81 ) 


/ i — 1 

The matrix is called a Jacobian matrix with a maximum size of 3 x 3. 

Thus Eq. (79) can be written as 

Since Eq. (79) represents a Taylor series approximation, the incremental 
correction Is used to obtain the next approximation 

AUj = AU]" 1 + U 1 (83) 

The relations in Eqs. (81) and (82) constitute the Newton-Raphson solution of 
Eq. (78). The iteration is continued until appropriate convergence criteria, 
discussed in the section Finite Element Formulation and Overall Scheme, are 
satisfied . 

f v 1-1 

The Jacobian matrix can be evaluated by finite difference 

perturbation techniques and placed in the following form 
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where 


M, 


[f m <AQ + dQ, AG, AJ) - f m (AQ, AG, AJ>] 
ml * dQ 

[f m <AQ, AG + dG, AJ) - f m <AQ, AG, AJ)] 




m2 


'm3 


dG 


[f m <AQ, AG, AJ + dJ> - f m (AQ, AG, AJ)] 


dJ 


m = 1,2,3 
dQ = (O.ODAQ, 
dG = (O.ODAG, 
dJ = (O.Ol)AJ. 

COUPLED CONTINUUM DAMAGE AND VISCOPLASTIC FORMULATION 


J 


( 84 ) 


(85) 


The nucleation of microcavities, and their growth and coalescence Into 
macroscopic cracks, Is generally the cause of material deterioration (material 
damage) such as decrease of strength, rigidity, toughness, stability, and 
residual life. Since the pionee - works of Kachanov (1958) and Rabotnov (1969), 
a new concept has been developed to Investigate the growth of microcavities and 
the mechanical behavior of damaged materials. This concept, called "continuum 
damage mechanics," represents the effects of distributed cavities In terms of 
certain mechanical variables. Since its notion hypothesizes that the effects 
of microcavities can be described by appropriate damage variables, such 
variables can be represented according to the same notion as that of stress, 
strain, or temperature field (Murakami, 1983). Therefore they are the same as 
the internal state variables in thermodynamical theories of constitutive 
equations . 

In this section we introduce a damage variable as an internal state 
variable and couple it with the constitutive and evolution equations of 
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Walker's viscoplastic theory. Our aim Is to demonstrate the potential 
advantage of the proposed Integration scheme when Investigating continuum 
damage models. Walker's damage model, which was selected for this exercise. 

Is first presented in differential form and then In uniformly valid asymptotic 
form. Note that a damaged state is not considered in the elastic constitutive 
equations . 

The concept of a damage variable was first proposed by Kachanov (1974) 
when he developed a mathematical model for evaluating creep rupture times. 
Cavity growth that results in a reduction of the net area Is assumed to be the 
principal mechanism of material damage. The damage state may be represented 
by an internal state variable <j/ such that y = 1 and >|> = 0 specify the 
undamaged initial state and the final rupture state, respectively. By taking 
the maximum effective stress a as the principal factor governing the 
progression of the damage, Kachanov formulated the evolution equation of the 
damage variable as follows: 



where A and r are material constants. Though Kachanov did not discuss the 
physical meaning of it may be interpreted as the ratio between the net 
area S n of a given section to that of the corresponding apparent area S 

S n 

* = f (87) 

The stress, which is magnified by the net area reduction, is called net stress 
and is defined 

a = — (88) 

n v 

As in the classical theories of creep, Eq. (88) can be generalized to 
mul tiaxial stress states. Assuming Isotropy of material and of material 
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damage, we, for demonstration purpose. Introduce t to Walker's viscoplastic 
model and derive the following equations: 
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( 86 ) 


(10) 


n], r\ 2 , n 3 , n 4 , ns, ng, m, A, and r are material constants. 

The UVAE form of Eqs. (89) and (90) is obtained in a fashion similar to 
that already described. However, the damage parameter is derived by directly 
integrating Eq. (86). The final form is as follows: 


a. j ( t + At ) = -^ Q 


§ v‘ * at> ♦ ( x * I “Kw 1 * 4t> 

* - f v t; - ( x * f p>ij'kk (t) ] 
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( 93 ) 
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Q^(t + At) = exp(-AG)|Q. j < t) - Q^J + n 2 Ae^— ^ J 
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K(t + At) = K 
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(95) 


i{/( t + At) 


_r, l/(r+l) 

1 - A(r + l)(t + At)a j 


(96) 


As mentioned earlier, the differential form of the unified viscoplastic 
formulation results in a mathematically stiff system of differential equations. 
When damage is incorporated, unstable behavior from the numerical integration 
tends to occur whether or not the explicit forward difference or the implicit 
backward difference method is used. This unstable phenomenon in the 
differential equations arises from the fact that the right side of Eq. (89) 
becomes very large and sensitive to the time-step increment as the damage 
parameter approaches zero. For the Integral form of this damage model, the 
factor y appears on the right side of Eq. (92), and this equation is the 
Intermediate term of Eq. (91) for the stress. When i|r approaches zero, AQ of 
Eq. (92) becomes large, and when AQ is substituted into Eq. (91), the stress 
decreases. There is no sign of numerical difficulty. Therefore, an unstable 
phenomenon should not be encountered if the proposed Integration scheme is used 
to integrate the continuum damage model. 

Both differential and UVAE forms of this coupled continuum damage and 
viscoplastic model have been Incorporated into a MARC finite element program. 
The results are presented in the section Numerical Examples and Comparisons. 
FINITE ELEMENT FORMULATION AND OVERALL SCHEME 

In the analysis of time-dependent constitutive relations, the formulation 
currently used is that for small strain theory; that is, material nonlinearity, 
only, is taken into consideration. In nonlinear finite element analysis, it is 
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most effective to use an Incremental formulation of the equation of motion. 

The global Incremental governing equtions are given as 

[KHAU} 1 = {AR} 1 (97) 

where {AU}' Is the vector of increments in nodal point displacements, {AR}^ 
is the vector of unbalanced force at iteration i, and [ K ] is the global 
stiffness matrix and is defined as 

[K] = f [8] T [C][B]dv (98) 

J v 

♦ 

where [ B ] is a strain-displacement transformation, and [C] is a stress-strain 
material property matrix. 

There are two well-known methods to represent the inelastic strain of 
constitutive relations governing each element at the local level and to 
assemble this information into the global level of Eq. (97). The first 
approach, the tangent stiffness method, combines elastic and inelastic strain 
characteristics at each increment directly into a tangent modulus [Cl, which 
is then supplied to the global equations and assembled into the global 
stiffness matrix [<]. This approach is commonly employed with rate-independent 
constitutive equations (i.e., plasticity). The second approach is called the 
initial strain method wherein the tangent modulus [C] is evaluated from the 
elastic material moduli only. The inelastic strain is carried to the global 
equations in the form of strain Increments {AeM- These strain increments are 
then assembled into a pseudo-load vector {AR*} which Is added to the right side 
of Eq. (97) and defined as 

{AR*} = [ [B] T [C] {Ae 1 }dv (99) 

J v 

This approach has been widely used with creep and unified viscoplastic 
constitutive equations and is employed throughout this work in the proposed 
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integration scheme presented in the previous sections. A flow chart describing 
the nonlinear finite element analysis with a global incremental iteration 
procedure is presented in Fig. 1 

At the local level the overall scheme with a Newtom-Raphson iteration is 
summarized in Fig. 2. The basic concept behind thfs scheme is using the UVAE 
equations of viscoplastic models, derived previously whi le ensuring that by 
taking only the first term of the expansion, the accuracy is obtained via 
Newton-Raphson iteration. Details are as follows: 

(1) With the Initial strain method, an inelastic strain Increment is 
assumed at the start of Iteration to be a devlatorfc strain Increment taken 
from the previous time step. However, for subsequent global iteration the 
deviatoric strain increment that is calculated fronrprevious global iteration 
is used. 

(2) The Initial guesses for A Q, AG, and AJ, or so-called local iteration 
vectors, are all assumed to be 0.1. These values acre judiciously chosen on 

the basis of experience. The values range between "0,1 and 3.0. For small and 
nonsevere loading problems, a high number is recommended; whereas for severe 
thermal and mechanical loading situations, a low number is more appropriate. 

The state variables oij, Qij, and K are then det erm ined by using Eqs. (71), 
(73), and (75). Whenever the local iteration vectoFls updated, these state 
variables must be recalculated. 

(3) To calculate the inelastic strain rate R,IJhe finite difference is 


employed by equating R to AR/At, where At is tM current time-step 

1/2 


increment and AR is determined as 


(2/3)Ae 


!. acM 

13 ijJ 


For the first 
i 


iteration, Ae^ is set to Ae^. In the subsequent.^ terations , Ac^ is set 
equal to Ae^ - ^AS.j/2p^). 
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(4) Once the state variables are known (based on the Initial guess of the 
local iteration vector), the functions in Eq. (77) can be evaluated to 
determine how good the guess is. These functions will be the right side vector 
of Eq. (82). Obviously, if the guess is good, these functions will be small, 
and the local iteration vector computed from Eq. (82) will also be small. This 
is the sign of convergence. 

(5) Next, the matrix is estimated from Eq. (84). The finite 

difference perturbation technique is based on 1 percent of the values of AQ, 
AG, and AJ. This proved to be a viable choice since no ill conditioning of 
of the matrix took place during iterations. 

(6) The iteration vector can now be corrected with Eq. (82) by inverting 

the matrix whose maximum size is only 3 x 3 as in Miller's model. For 

elasto-pl asti c-creep model, the matrix reduces to a scalar. Inverting 

these matrices costs nothing; this is a tremendous advantage when analyzing 
large-scale problems. Once the iteration vector is corrected, the new values 
of AQ, AG, and AJ, as well as the state variables , Q.^, K, e^, and 
ejj, are subsequently updated. 

(7) One of the most important parts of this iterative scheme is the 
convergence criteria. In order for the algorithm to be effective, realistic 
criteria should be utilized for terminating the iteration process. At the end 
of each iteration, the solution that has been obtained should be checked to 
see whether it has converged within preset tolerances or whether the iteration 
is diverging. If the convergence tolerances are too loose, inaccurate results 
are obtained; if the tolerances are too tight, excessive computational effort 
is wasted for needless accuracy. Three convergence criteria are incorporated 
into the proposed integration scheme. First, the Iteration vector convergence 
criterion is defined as 
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< ETOLB 



where ETOLB is a convergence tolerance set equal to 0.01. The term Uj * 2 

(k) 

Is the two-norm of Iteration vector In the first iteration. The 2 

the two-norm of correction vector in the subsequent iterations. If the 
criterion is satisfied, the solution is obtained. There is no need to check 
other criteria. However, if it is not satisfied, two Cauchy stress convergence 
criteria need to be satisfied for the solution to converge. The first Cauchy 


stress convergence criterion is defined as 


(k+1 ) (k) 

Oi i - <7; i 


ij II 2 


(k-1 ) 
a ij 


This is a criterion to prevent any unnecessary iterations. If it is satisfied, 
the second Cauchy sress convergence criterion is checked; it is defined as 


follows: 

< CTOL 

where CTOL is set equal to 0.005. This Is a fairly tight tolerance. 

The Cauchy stress was the only state variable selected for convergence 
checks because Cauchy stress is the only state variable needed at the global 
level. In contrast, the back and drag stresses have never been used at the 



global level . 

NUMERICAL EXAMPLES AND COMPARISONS 

To demonstrate the numerical behavior of the new integration scheme, it 
was coded into subroutine HYPELA (see Appendix D), which is written in FORTRAN 
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and is interfaced with the MARC finite element program. The integral form of 
Walker's, Krieg, Swearengen, and Rohde's (KSR's) and Miller's constitutive 
models, as presented herein were incorporated into one subroutine. The 
subroutine is written in a very efficient and versatile way, as all three 
models are included and tied to one integration scheme. A HYPELA subroutine 
containing the differential form of Walker's model with a SAFE integration 
scheme was taken from Cassenti (1983b) and used primarily for comparison with 
the proposed UVAE integration scheme. Similar subroutines containing the 
differential form of KSR's and Miller's models with a SAFE scheme were also 
written and used in the same fashion; these are not provided in this report. 

For Walker's damage model presented in the section Coupled Continuum Damage 
and Viscoplastic Formulation, minor modifications that are needed can be made 
with relative ease to subroutine HYPELA, for both differential and integral 
forms. Hence, it is not reproduced in this work. All the analyses were 
performed on the Cray-XMP super computer at NASA Lewis Research Center. 
Comparions are based on the number of seconds of Central Processing Unit (CPU) 
time used. 

Hysteresis Loop for Hastelloy-x Under Thermomechanical Loading at 1600 °F 

An axi symmetric finite element model was used to simulate one quarter of 

a solid specimen made of Hastelloy-x metal, which is being used for jet engine 

combustor liners. The material constants of Hastelloy-x at 1600 °F for 

Walker's, KSR's, and Miller's models are given in Cassenti (1983a). The cyclic 

• 3 1 

response is governed by these parameters: strain rate e = 3.87x10 sec 

and strain limit of 0.006 in. /in. One full cycle of load, consisting of three 
loading portions, was imposed as follows: portion 1 - loading, 0 < e < 0.006; 
portion 2 - unloading, 0.006 > e > -0.006; and portion 3 - loading, 

-0.006 < e < 0.006. For each model, three different total time steps, with 
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equal step sizes, were specified (i.e., 20, 40, and 80 time steps). Our 
purpose was to study the stability and accuracy of the algorithm as time-step 
size was increased. 

Figure 3 shows the hysteresis loop of Walker's model for 80 time steps. 
Three different runs are plotted: SAFE integration with ERR0R1 = 0.0001; 

SAFE integration with ERROR 1 = 0.00001; and proposed UVAE integration with 
ETOLB = 0.01 and CT0L = 0.005. The results are almost identical for all three 
runs. For 40 and 20 time steps, similar results were obtained, as shown in 
Figs. 4 and 5, respectively. However, comparison of CPU time (see Table I) 
indicates that the proposed UVAE scheme is more efficient computationally as 
the step size increases. Of course, greater accuracy is attained because of 
the iterative nature of the scheme, as compared to the SAFE scheme which is a 
noniterative type. 

Figure 6 shows the results of the same three runs with KSR's model for 
40 time steps. Results are again nearly identical. By comparing the CPU 
times summarized in Table I, conclusions similar to those for Walker's model 
can be drawn. For Miller's model, the amounts of CPU time are quite different. 
Table I shows that the proposed UVAE scheme consumes 2.5 times more CPU time 
than does the SAFE scheme with ERR0R1 = 0.0001 for the case of 80 time steps. 
However, as the number of time steps decreases or the size of time steps 
increases, the efficiency becomes comparable. Good accuracy is obtained for 
both schemes as shown in Fig. 7. Because this model has a very stiff region, 
the user has a tendency to be more conservative and specifies a tight tolerance 
for the SAFE scheme with ERR0R1 = 0.00001; then CPU time is three to four times 
higher than with ERR0R1 = 0.0001, as can be seen in Table I. The noniterative 
nature of the scheme lures the user to be conservative; however, this does not 
happen with the UVAE scheme since accuracy is always assured through iteration. 
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Thermomechanical Fatigue Loops 

Thermomechanical fatigue loops (TMF) are typical loading histories 
experienced by Hastelloy-x material in jet engine combustor liners at elevated 
temperatures. Linder such conditions, both mechanical load (in the form of 
imposed strain) and temperature are subjected to large changes as a function 
of time. In order to predict the life of combustor liners realistically, the 
analyst must have a precise knowledge of the stress-strain hysteresis behavior 
at the critical fatigue locations corresponding to the aforementioned loadings. 
The purpose of considering the TMF is twofold: the first is to demonstrate the 

capability of the proposed UVAE scheme in handling noni sothermal loadings, and 
the second is to assess the predictive capability of the Walker and KSR models, 
based on the proposed scheme, as compared to the experimental data reported in 
Cassenti (1983a, 1983b). 

Considered herein is the case of an open nonsymmetri cal TMF cycle as 
shown in Fig. 8. The temperature varies sinusoidally from 950 to 1750 °F, 
with a temperature hold at 1750 °F for 40 sec; the strain, which also varies 
sinusoidally, holds -0.43 percent for the same period. The total number of 
time steps used for all analyses was 56. 

The results of using Walker's model for a SAFE scheme and using the UVAE 
integration scheme are presented in Figs. 9 and 10 along with the experimental 
results. Notice that the proposed scheme gives better results than the SAFE 
scheme (especially during steady-state conditions) when both are compared to 
the experimental results. The superiority of the proposed scheme's results 
are even more obvious in the KSR model results shown in Figs. 11 and 12. 
Comparisons of CPU time for both schemes and both models are shown in Table II. 
The new scheme utilizes only 5 percent more CPU time for Walker's model, and 
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25 percent more for KSR's model, than does the SAFE scheme. However the new 
scheme's accuracy is undeniably better. 

Annular Combustor Liner Test Rig 

In this example, a large-scale analysis of a combustor liner was used to 
demonstrate the effectiveness of the UVAE scheme. The combustor liner is a 
cylindrical part of a gas turbine engine that was radiantly heated in the 
structural component response rig in the test liners. A photograph of the 
conventional test liner is shown in Fig. 13(a). The test liner, of sheet- 
metal seam-welded louver construction, is a nickel-base superalloy material, 
Hastelloy-x. The eight louvers are segments of an outer annulus of a 
combustor liner. The test liner has an inside diameter of approximately 
50.8 cm (20 in.). Circumferential arrays of cooling holes cool the louver 
lips. Louvers 4 to 6 (see Fig. 13(b)) are the active test louvers, that is, 
the location where the heat flux to the test liner is considered to be 
relatively flat. 

A typical engine's mission cycle (takeoff, cruise, landing, and taxi) of 
3 to 4 hr was simulated in 2.2 min. This thermal cycle time is broken up into 
four segments (see Fig. 14): a 6-sec ramp up from minimum to maximum power; a 

1-min hold time at maximum power; a 6-sec ramp down from maximum to minimum 
power; and a 1-min hold time at minimum power. Cyclic surface temperatures at 
two potentially critical failure locations (the seam weld and the knuckle) on 
the liner of louver 5 are plotted in Fig. 15. These data were used in the 
heat transfer analysis, with MARC code, as boundary conditions (thermal loads). 
Details of this analysis can be found in Thompson and Tong (1986). 

The output of the heat transfer analysis was used as input to the 
structural analysis program. A three-dimensional solid finite element model 
of louver 5, consisting of 546 elements and 1274 nodes, is shown in Fig. 16. 
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Appropriate boundary conditions are assumed. Walker's model was used to 
perform the analysis for both the SAFE and the UVAE integration schemes. 
Because of the large amount of CPU time consumed, the results of only 10 time 
steps are presented. In Figs. 17 and 18, hoop stress and strain at the seam 
weld and at the knuckle, respectively, are plotted. The results are in good 
agreement. Discrepancies are within 3 percent at the knuckle and 10 percent 
at the seam weld. However, a question arises about whether the SAFE scheme is 
accurate for a large-scale analysis with complex loading histories, since it 
is a subi ncremental noniterative approach. The error that occurs in each time 
step of a large-scale analysis may be sizable and cumulative. Because of the 
iterative approach of the proposed UVAE scheme, error is not accumulated. 
Comparison of CPU times shows that the UVAE scheme (271 sec) has an 8-percent 
advantage over the SAFE scheme (292 sec) when only 10 time-step increments are 
used . 

Continuum Damage Behavior During Creep Rupture Test 

In the section Coupled Continuum Damage and Viscoplastic Formulation, a 
damage model was incorporated into Walker's viscoplastic model in both 
differential and Integral forms. The damage parameter ifr was introduced. To 
test these models numerically, subroutine HYPELA was slightly modified. The 
finite element model is the same as for the Hastelloy-x hysteresis loop. It 
is first loaded to stress, which saturates at a value of 7500 psi throughout 
the analysis. The values of A and r were chosen as 6 . 208 19891x1 0 — 26 anc j 
5.4, respectively, to provide verification of the numerical scheme. The damage 
parameter y was initially set equal to 1 and diminished toward zero, as shown 
in Fig. 19. No numerical difficulty was encountered. However, for the SAFE 
scheme, a breakdown occurred at y = 0.58, even though a very small time step 
was specified. At creep rupture the strain was 0.0058, and the rupture time 
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was 3300 sec (see Fig. 20). This example demonstrated that the proposed UVAE 
integration scheme possesses a tremendous advantage In the analysis of 
continuum damage mechanics. 

SUMMARY OF RESULTS 

A new uniformly valid asymptotic Implicit integration algorithm for 
elasto-plastlc-creep and unified viscoplastic theories, Including continuum 
damage, is proposed and demonstrated through a user subroutine of the MARC 
commercial finite element code. Based on the results obtained, the following 
characteristics of the proposed algorithm can be stated: 

1. The algorithm is iterative without a high computational cost. 

2. The algorithm is stable for large time Increments. 

3. The results obtained are less user-dependent. 

4. The algorithm is simple, easy to implement, and well suited for finite 
element appl 1 cations . 

5. Under complex loading histories, including multiaxial behaviors, the 
algorithm is accurate and efficient. 

6. The algorithm was shown to possess a tremendous advantage in continuum 
damage mechanics. 

7. The algorithm is suitable for large scale multiaxial problems. 
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APPENDIX A 


WALKER'S THEORY 


(1) Differential Form 
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and X, jj, Q.j, n, m, n^, n 3> n g , n ? , K ] , and < 2 are material constants and 
depend on temperature. 
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APPENDIX B 


KRIEG, SWEARENGEN, AND ROHDE'S (KSR) THEORY 
(1) Differential Form 


•1 _ / 

c,j "V 
1 


) [f(! s n - a n)(l s n - a n)l 

K 


1 / 2 ' 


(I !u - 


u hi 


j 


[id s ij - «h)(! 


1J/Y2 1j 1J 




, 1/2 


* 1/2 

Q ij = A 1 S iJ - A 2 Q..(! Q pq Q pq ) [exp(A 3 § Q pq Q pq ) - l] 


K = K, 


(2) Uniformly Valid Asymptotic Expansion Form 

, 1j<t ♦ it) - f a,j<t * it) , (x ♦ f u)s ljCkk (t * it) 

i<t) - § Q 1 J (t) - (x + f l i)6 ij e kk (t)]exp(-AQ) 




( B1 ) 

( B2 ) 
(B3) 


where 


( 2 P -f Mn Ae,,,, - \ Afl 


\H - exp(-AQ)! 

13 3 ^ u ij uc -kk - 3 u “ij;L aq ■ J 


Q 1 j(t + At) = exp(-AG)0 1 j(t) + A, Aejj[— ~ 


K = K r 


(B4) 

(B5) 

(B6) 


3p At *,*. . j.v l-< 1 /n) 

AQ - koVao R(t + At) 

1 /2 

AG = !\J\ f) Q ) fexp^A. | Q Q ) - 1 1 At 

2Y3 pq pq / L K V 3 3 pq pq/ J 

and X, p, n, A 1 , A 2 , A^, and K Q are material constants which depend on 
temperature . 


( B7 ) 


( B8 ) 
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APPENDIX C 


MILLER'S THEORY 


<1) Differential Form 


e.j = B0'\s1nh 


3(2 S i.i Q i.l)(I S i.i 

- 

K 



,3/2^1 n 


(iiu-ii) 


If s .j - s lj - n u)] 


1/2 


Q.j = - H ] B6’ <sinh 


‘lG Q ij fi 1j) 


1/2 


V\n 


11 


(f Vu) 


1/2 


K = H 2 R 


1/2 A 

C 2 * (I Q 1 3 Q 1j) -A7< 3 - H 2 C 2 BO'[ S inh(A 2 K 3 )] 


where 


6 1 = exp(- 21) for T > 0.6 T ( 
/0. 6 T, 


m 


0' = exp 


(- -V—) 

V 0.6 kT ) 


2,n 


m 


for T < 0.6 T, 


m 


(2) Uniformly Valid Asymptotic Expansion Form 

ij<t * «> - § Q,j<t t at) ♦ (x * | ii)«,j« kk (t * At) 

+ [o. j<t) - | fl^(t) - (x + | p)s ijE|ck (t)]exp(-AQ) 

♦(* ^13 - f ^ S U - 1 ^ia)[^gr^] 

Q i j(t + At) = exp(-AG)Q..(t) + H 1 Ae]j P ~ e *g ( ~ A -~] 

exp(-AJ) + H 2 C 2 Ar[ ] ~ e ^ ~ A — ] 


K(t + At) = K q + K(t) - K q j 


(Cl ) 

(C2) 

(C3) 

(C4) 

(C5) 


(C6) 

(C7) 

(C8) 



where 



At 



2/3 



and T m> n, H 2> A 1 , A 2 , B, C 2 , Q*, and k are material constants which 
Independent of temperature. The material constants X, \i, H 1 , Kg, and 
depend on temperature; T Is the temperature in Kelvin. 


(C9> 


(CIO) 


(Cl 1 ) 

are 

9 ' 
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APPENDIX D 


SUBROUTINE HYPELA 

SUBROUTINE HYPELA (D , G , E , DE , S , TEMP , DTEMP , NGENS , N , NN , KC , MAT , NDI , 

1 NS HEAR) 

C 

C A NEW VALID ASYMTOTIC INTEGRATION SCHEME FOR 3 VISCOPLASTIC MODELS 
C 

C MOD = 1, WALKER’S MODEL 

C MOD = 2, KSR’S MODEL 

C MOD = 3, MILLER’S MODEL 

C 

C*K*ttt**»n THIS SCHEME IS WRITTEN BY * A. CHULYA "«*******«*»***.** 
C ALL RIGHTS RESERVED 

C 

DIMENSION D (NGENS , NGENS) , G (NGENS) , E (NGENS) , DE (NGENS) , S (NGENS) 
DIMENSION TEMP ( 1 ) , DTEMP ( 1 ) 

DIMENSION SIGB(6) ,0MEGB(6) ,CB(6) ,SIGE(6) , OMEGE(6) ,CE(6) 

DIMENSION DC (6) , DET (6) , OMEGI (6) 

DIMENSION DSIGIN(6) ,DS(6) ,AB(6) 

DIMENSION F (2 , 3) ,BUP(3) ,DCTEMP(6) ,TlSIG(6) ,FM(3,4) 

COMMON/AKEV/KEVIN 

COMMON/FAR/DUM, INC 

COMMON/CDC/DUMMY(I8) , NCYCLE 

* * f * r 

..... SECOND INVARIANT FUNCTION 

SINV(A,B,C,D,E,F)=(A*A+B*B*C«C+2. » (D*D+E»E+F*F) ) *2 . /3. 

C 

C USERS SELECT THE VISCOPLASTIC MODEL 
MOD = 1 

C 

IF (MOD. LE. 3) GO TO 9 
WRITE (6,4711) 

4711 FORMAT (’ MODEL SELECTED IS INVALID - SOLUTION STOP ’) 

STOP 

C. ....DETERMINE IF PLANE STRESS, PLANE STRAIN , AXISYMMETRIC , OR 3-D 
C.»*.*KELTYP=1 FOR PLANE STRAIN AND AXISYMMETRIC PROBLEMS 
C****»KELTYP=2 FOR PLANE STRESS PROBLEM 
C*....KELTYP=3 FOR 3-D PROBLEM 

9 IF(NDI.EQ.3.AND.NSHEAR.EQ.l) KELTYP=1 
IF (NDI . EQ . 2 . AND . NSHEAR . EQ . 1) KELTYP=2 
IF (NDI . EQ . 3 . AND . NSHEAR . EQ . 3) KELTYP=3 
C.....SET UP CONSTANTS 
MAXIT=25 
NELPR=1 
IPR=1 
NPRIN= 1 
SFTEMP-936.2 
C. . . .SET IT TOLERANCE 
ETOLB 0.01 
CTOL = 0.005 

PUT STRESSES AT BEGINNING OF MARC INCREMENT INTO SIGB ARRAY ACCORD 
C..*..TO ELEMENT TYPE 

GO TO (801 , 802 , 803) , KELTYP 

801 CONTINUE 
SICB(1)=S(1) 

SIGB(2)=S(2) 

SIGB(3)=S(3) 

SIGB(4)=S(4) 

SIGB (5) =0. 

SIGB (6) =0. 

GO TO 900 

802 CONTINUE 
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SIGB(1)=S(1) 

SIGB(2)=S(2) 

SIGB(3)=0. 

SIGB(4)=S(3) 

SIGB(5)=0. 

SIGB(6)=0. 

GO TO 900 

803 DO 804 J=1 , 6 
SIGB(J)=S(J) 

804 CONTINUE 
900 CONTINUE 

C*****INITIALIZE STATE VARIABLES ON FIRST ENTRY TO SUBROUTINE. ON SECOND 
C AND SUBSEQUENT ENTRIES SKIP INITIALIZATION. 

KEVIN=INC+NCYCLE 
IF (KEVIN. NE.O) GO TO 3 
TEMP ( 1 ) = SFTEMP 
DO 2 J=2, 15 
TEMP ( J) =0. 

2 CONTINUE 

3 CONTINUE 

C***SET STARTING VALUES OF STATE VARIABLES DURING PRESENT MARC INCREMENT 
DEG=TEMP(1) 

TB=TEMP(2) 

RB=TEMP (3) 

IF (MOD. EQ. 3) AKB=TEMP(16) 

DO 104 KA=1 , 6 
J = KA+3 

OMEGB (KA) =TEMP ( J) 

CB(KA)=TEMP(J+6) 

104 CONTINUE 

C SET TEMPERATURE AND TIME SUBINCREMENTS 

DDEG=DTEMP ( 1 ) 

DT=DTEMP(2) 

C PUT SUBINCREMENTS OF TOTAL STRAIN INTO ARRAY DET ACCORDING 

C TO ELEMENT TYPE 

C 

GO TO (61,62,63) , KELTYP 

61 CONTINUE 
DET(l) = DE(1) 

DET (2) = DE(2) 

DET (3) = DE(3) 

DET (4) = 0. 5*DE(4) 

DET (5) = 0. 

DET(6) = 0. 

GO TO 71 

62 DET ( 1 ) = DE(1) 

DET (2) = DE (2) 

DET (3) = -DET(l) -DET (2) 

DET (4) = 0.5*DE(3) 

DET (6) = 0. 

DET(6) = 0. 

GO TO 71 

63 CONTINUE 

DO 64 J=1 , 6 
FAC=1 . 

IF(J.GT.3)FAC=0.5 
DET(J) = FAC*DE(J) 

64 CONTINUE 
71 CONTINUE 

C****+SET INITIAL GUESS FOR EQUILIBRIUM STRESS AT END 
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C***+*0F SUBINCREMENT EQUAL TO EQUILIBRIUM STRESS AT 
(^♦♦♦♦♦BEGINNING OF MARC INCREMENT 
C 

DO 2000 J=l , 6 
OMEGE(J) = OMEGB(J) 

2000 CONTINUE 

♦♦♦ ♦♦ASSUME INITIAL GUESS FOR INELASTIC STRAIN IN FIRST ITERATION 
♦♦♦♦♦EQUAL TO DEVIATORIC STRAIN INCREMENT 

DV0L=DET(1)+DET (2) +DET (3) 

DO 72 J=1 , 6 
ALPHA = 1. 

IF(J.GT.3) ALPHA-0 

DC(J) = DET(J) - ALPHA *DV0L/3. 

72 CONTINUE 

♦♦♦♦♦COMPUTE INELASTIC STRAINS AT END OF FIRST SUBINCREMENT 


DO 7125 J=1 , 6 
CE(J) = CB(J) *DC(J) 

7125 CONTINUE 

.START INTEGRATION 

♦♦♦♦♦COMPUTE TEMPERATURE DEPENDENT MATERIAL CONSTANTS 
DEGM = DEG + 0.5*DDEG 
GO TO (41,42,43) , MOD 

41 CALL C0NST1 (DEGM,EE,ANU, AK0,ANIN,AM,AN1 , AN2 , AN3 , AN4 , 

1 AN5 AN6 , AN7 , OMEGZ , AN , AL AM , AMU , Cl , C2 , C3 , C4 , C5) 

C**+**SET INITIAL VALUES OF EQUILIBRIUM STRESS 

DEN0M=SINV(CE(1) ,CE(2) ,CE(3) ,CE(4) ,CE(5) ,CE(6)) 

AB(l)= D 0MEGz!2 E »0MEGZ^(CE(l)*CE(l)+CE(4)»CE(4)-CE(6) ‘CE(6 'i 

*AB(2)=-0MEGZ+2. ♦OMEGZ* (GE (4) +CE(4) *CE(2) »CE(2) +CE(5) ♦CE(5) 

1 AB(3)=-0MEGZ+2. ♦OMEGZ* (CE(6) +CE(6) *CE(5) ♦CE(5) *CE(3) *CE(3) * 

AB(4)=2 / .0MEGZ* (CE(1) *CE(4) *CE(2) ♦CE(4) +CE(5) ♦CE(6) + 1 .E-30)/ 
1 DFNOM 

AB(5)=2 . ♦OMEGZ* (CE(4) *CE(6) *CE(2) *CE(5) +CE(3) .CE(5) +1 . E-30) / 

1 AB (6^=2 . *OMEGZ* (CE ( 1 ) *CE (6) +CE (4) +CE (5) *CE (3) *CE (6) * 1 . E-30) / 
1DEN0M 

ABSUM=AB(1)*AB(2)*AB(3) 

DO 7134 .1=1,6 
ALPHA=1 . 

IF (J. GT . 3)ALPHA=0. 

OMEGI ( J) =AB ( J) - ALPHA* ABSUM/3 . 

7134 CONTINUE. 

GO TO 69 

42 CALL C0NST2 (DEGM , EE , ANU , AKO , ANIN , A 1 , A2 , A3 , A 4 , 

1 A5,AN,ALAM,AMU,C1,C2,C3,C4,C5) 

C* * * * *SET INITIAL VALUES OF EQUILIBRIUM STRESS 
DO 7135 J=1 ,6 
OMEGI ( J) =0 . 

7135 CONTINUE 
GO TO 69 

43 CALL C0NST3 (DEGM , EE , ANU , AKO , AN , A1 , A2 , HI , H2 , Z2 , BP , 

1 ALAM, AMU,C1 , C2 ,C3 ,C4 ,C5) 
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DO 7136 J=1 ,6 
OMEGI (J)=0. 

7136 CONTINUE 
C 

IF (KEVIN. EQ.O) AKB=AKO 
C 

C COMPUTE INITIAL DR AND SAVE IN DRC 

C 

69 DR = SINV (DC (1), DC (2), DC (3), DC (4), DC (5), DC (6)) 
DR = SQRT(DR) 

IF(DR.LE.l.E-lO) DR=1.E-10 
NIT = 0 
C 

C ASSUME INITIAL GUESSES 

C 

DQ = 0.1 
DG = 0.1 

IF (MOD. EQ. 3) DJ = 0.1 
C 


C START ITERATION LOOP 

C 

73 CONTINUE 
NIT = NIT+1 

C CALCULATE THE BACK STRESS 

IF(MOD.EQ.l) DCON = AN2 
IF (MOD. EQ. 2) DCON = A1 


IF (MOD. EQ. 3) DCON = HI 

CALL OMEGAR (DG , DCON , OMEGI , OMEGB , OMEGE , DC) 

C 

C CALCULATE THE DRAG STRESS 

CALL KAPPAR (D J , H2 , Z2 , AKO , AKB , AKE , DR , MOD) 

C 

C CALCULATE STRESS 

CALL SIGMAR (DQ , KELTYP , DET , AMU , ALAM , DC , OMEGB , OMEGE , 

1 SIGB,SIGE,DVOL) 

STNORM=SQRT (SIGE (1) * *2+SIGE (2) **2+SIGE (3) * *2+SIGE (4) * *2 
1 +SIGE(5) **2+SIGE(6) **2) 

C 

DO 435 K=1 ,6 
ALPHA = 1. 

IF (K . GT . 3) ALPHA=0 . 

DCTEMP(K) = DC(K) 

DC (K) = (ALPHA*ALAM*DV0L+2 . *AMU*DET (K) -SIGE (K) +SIGB (K) ) / 
ft (2. * AMU) 

435 CONTINUE 
DRTEMP = DR 
C 

C COMPUTE DELTA R FOR NIT >= 2 

C 

IF (NIT.EQ. 1) GO TO 444 

CALL DELR (NIT , DVOL , ALAM , AMU, DET , SIGB , SIGE , 

1 OMEGE, AKE, AN, DT, DR) 

444 RDOT = DR/DT 

GO TO (541,542,543), MOD 

541 CALL EVALF1 (RDOT,3,AN,AM,AN3,AN6,AMU,DT,F,OMEGE, 

1 AKE,AKO,DQ,DG) 

GO TO 544 

542 CALL EVALF2 (RD0T,3,AN,A2,A3,AMU,DT,F,0MEGE, 

1 AKE, AKO, DQ, DC) 

GO TO 544 
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543 CALL EVALF3 (RDOT , 4 , AN , BP , Hi , H2 , A1 , A2 , Z2 , AMU , DT , FM, OMEGE , 

1 AKE,AKO,DQ,DG,DJ) 

544 IF (NIT. EQ. 1) GO TO 405 
C 

C CONVERGENCE CHECK 

IF (BNORM . LE . (ETOLB+DNORM) ) GO TO 909 
C 

SUM1=0.0 
DO 2002 1=1,6 

SUM1=SUM1+ (SIGE(I) -TlSIG(I ) )**2 
2002 CONTINUE 

SUM1=SQRT(SUM1) 

IF (NIT. LE. 2) GO TO 404 
C 

IF(SUM1 . GT.SUM2) GO TO 404 
IF (SUM2 , GT . CTOL*T2NORM) GO TO 404 
C 

C. . , . .UPDATE INELASTIC STRAIN C AT T+DT 
C 

909 DO 911 1=1,6 

CE(I) = CB(I) + DC (I) 

911 CONTINUE 
C 

C COMPUTE INELASTIC STRESS INCREMENT ACCORDING TO ELEMENT TYPE 

C 

GO TO (809,810,809) ,KELTYP 

809 DO 812 J=1 , 6 
DSIGIN(J)=-2 . *AMU*DC (J) 

812 CONTINUE 
GO TO 902 

810 DO 813 J=l,6 
ALPHA = 1. 

IF(J.GT.3) ALPHA=0. 

DSIGIN ( J) =ALPHA*2 . *AMU*ALAM+DC (3) / (ALAM+2 . * AMU) -2 . ♦ AMU*DC ( J) 

813 CONTINUE 
902 CONTINUE 

GO TO 420 
C 

404 IF (NIT . GE . MAXIT) GO TO 1991 
C 

405 DO 2001 1=1,6 
TlSIG(I) = SIGE(I) 

2001 CONTINUE 

T2N0RM = T1N0RM 
T1N0RM = STNORM 
SUM2 = SUM1 
C 

K = 2 

IF(yOD.EQ.3) K = 3 

DO 399 J=1,K 

GO TO (421,422,423) , J 

421 TEMPD = .01*DQ 

IF (ABS (TEMPD) . LE . 1 . E-8) TEWPD = l.E-4 
DQ1 = DQ + TEMPD 
DG1 = DG 

IF (MOD. EQ. 3) DJI = DJ 
GO TO 430 

422 TEMPD = .01*DG 

IF (ABS (TEMPO) .LE.l. E-8) TEMPD = l.E-4 
DQ1 = DQ 
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DG1 = DG + TEMPD 
GO TO 430 

423 TEMPD = ,01*DJ 

IF (ABS (TEMPD) .LE.l.E-8) TEMPD = l.E-4 
DG1 = DG 

DJI = DJ + TEMPD 
C 

C CALCULATE THE BACK STRESS 

IF(MOD.EQ.l) DCON = AN2 
IF (MOD. EQ. 2) DCON = A1 
IF (MOD. EQ. 3) DCON = HI 

430 CALL OMEGAR (DG1 , DCON, OMEGI , OMEGB , OMEGE , DCTEMP) 

C 

C CALCULATE THE DRAG STRESS 

IF (MOD. EQ. 3) CALL KAPPAR (D Jl , H2 , Z2 , AKO , AKB , AKE , DRTEMP , MOD) 
C 

C CALCULATE STRESS 

CALL SIGMAR (DQ 1 , KELTYP , DET , AMU , ALAM , DCTEMP , OMEGB , OMEGE , 

1 SIGB,SIGE,DVOL) 

.COMPUTE DELTA R FOR NIT=2 AND 3 

IF (NIT.EQ. 1) GO TO 560 
CALL DELR (NIT , DVOL , ALAM , AMU, DET, SIGB , SIGE , 

1 OMEGE, AKE, AN, DT,DR1) 

.COMPUTE RATE OF R 

RDOT = DR1/DT 
GO TO 561 

560 RDOT = DRTEMP /DT 

561 GO TO (571,572,573), MOD 

571 CALL EVALF1 (RDOT , J , AN , AM , AN3 , AN6 , AMU , DT , F , OMEGE , 

1 AKE,AK0,DQ1,DG1) 

GO TO 574 

572 CALL EVALF2 (RDOT , J , AN , A2 , A3 , AMU , DT , F , OMEGE , 

1 AKE,AK0,DQ1,DG1) 

GO TO 574 

573 CALL EVALF3(RD0T,J,AN,BP,H1,H2,A1,A2,Z2,AMU,DT,FM, OMEGE, 

1 AKE, AK0,DQ1 ,DG1 ,DJ1) 

C 

574 IF (MOD. EQ. 3) GO TO 580 

F(1 , J) = (F(l , J) -F(l ,3))/TEMPD 
F(2, J) = (F (2 , J) -F (2,3)) /TEMPD 
GO TO 399 

580 FM(1 , J) = (FM(1 , J) -FM ( 1 , 4) ) /TEMPD 
FM(2, J) = (FM (2 , J) -FM (2 , 4) ) /TEMPD 
FM(3, J) = (FM (3 , J) -FM (3 , 4) ) /TEMPD 
C 

399 CONTINUE 
C 

IF (MOD. LE. 2) CALL INVER2 (F , BUP) 

IF (MOD. EQ. 3) CALL INVER3(FM,BUP) 

C 

BNORM=SQRT (BUP (1) *BUP (1) +BUP (2) *BUP (2) +BUP (3) *BUP (3) ) 

IF(NIT.NE.l) GO TO 469 

DNORM=BNORM 

IF (MOD . LE . 2) ZNORM = SQRT(F(1 ,3) *F(1 ,3) +F(2,3) *F(2, 3)) 

IF (MOD . EQ . 3) ZNORM = 

1 SQRT (FM (1,4) *FM (1,4) +FM (2,4) *FM (2 , 4) +FM (3 , 4) *FM (3,4)) 
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c 

C UPDATE DQ, DG ft DJ 

C 

469 DQ = DQ + BUP(l) 

DG = DG + BUP(2) 

IF (MOD. EQ. 3) DJ = DJ + BUP(3) 

C 

GO TO 73 

C*****END OF ITERATION LOOP 

C*****PUT ELASTICITY MATRIX IN D AND INELASTIC STRESS INCREMENT IN G 
420 GO T0(814,815,816) .KELTYP 

814 CONTINUE 

DO 817 J=l,4 
DO 817 K=l,4 
D(J,K)=0. 

817 CONTINUE 

DO 818 J=1 , 3 
DO 818 K=l,3 
ALPHA=0 . 

IF(J.EQ.K) ALPHA=1 . 

D ( J , K) =C5+ALPHA*C3 

818 CONTINUE 
D(4,4)=C4 
GO TO 903 

815 CONTINUE 
D (1 , 1)=C2 
D(1 ,2)=C1 
D(2, 1)=C1 
D ( 1 ,3) =0 . 

D(3, 1)=0. 

D(2,2)=C2 

D(2,3)=0. 

D(3,2)=0. 

D(3,3)=C4 
GO TO 903 

816 CONTINUE 

DO 819 J=1 ,6 
DO 819 K=1 , 6 
D(J,K) =0. 

819 CONTINUE 

DO 820 J=1 ,3 
DO 820 K=l,3 
ALPHA=0 . 

IF(J.EQ.K) ALPHA=1 . 

D ( J, K) =C5+ALPHA*C3 

820 CONTINUE 
D(4,4)=C4 
D(5,5)=C4 
D(6,6)=C4 

903 CONTINUE 

DO 821 J=1 ,NGENS 
G(J)=DSIGIN(J) 

821 CONTINUE 

C+****COMPUTE STRESS AT END OF MARC INCREMENT 
DO 822 J=1 , NGENS 
SUM=0. 

DO 823 K=1 , NGENS 
SUM=SUM+D(J,K) *DE (K) 

823 CONTINUE 

DS(J)=SUM+G(J) 
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822 CONTINUE 

C*****PUT STATE VARIABLE INCREMENTS IN TEMP ARRAY FOR NEXT MARC INCREMEN 
DTEMP (3) -DR 

IF (MOD. LE. 2) GO TO 850 
DTEMP (16) =AKE-AKB 
IF (KEVIN. EQ.O) TEMP(16)=AK0 
850 DO 923 KA=1,6 
J=KA+3 

DTEMP (J) =OMEGE (KA) -TEMP ( J) 

DTEMP ( J + 6) =CE (KA) -TEMP ( J+6) 

923 CONTINUE 

IF (IPR.EQ.O) GO TO 12 
IF (NELPR . NE . N) GO TO 12 
IF (NN.NE.NPRIN) GO TO 12 
IF (NCYCLE.EQ.O) NWALK=0 
NWALK = NWALK + 1 
NQ = NWALK-2*NCYCLE 
NQQ=NCYCLE-1 
WRITE (6 ,20) INC 
20 FORMAT ( ’ INCREMENT’ ,15) 

WRITE (6, 750) NIT 
750 FORMAT (’ ITERATIONS’ ,15) 

WRITE (6, 753) N,NN 

753 FORMAT (’ ELEMENT ’, 15 , ’ INTEGRATION POINT’, 15) 

IF (NQ. EQ.O) WRITE (6, 23) NQQ 
IF(NQ.GT.O) WRITE (6, 39) NCYCLE 

23 FORMAT (55H VALUES OF PARAMETERS DURING SOLUTION OF RECYCLE NUMBER, 
115) 

39 FORMAT (55H VALUES OF PARAMETERS DURING ASSEMBLY OF RECYCLE NUMBER, 
115) 

WRITE (6, 29) 

29 F0RMAT(18H STRAIN INCREMENTS) 

WRITE (6,30) (DE (J) , J=1 , NGENS) 

30 F0RMAT(1P6E15.6) 

WRITE(6,31) 

31 FORMAT (18H STRESS INCREMENTS) 

WRITE (6 , 30) (DS ( J) , J=1 , NGENS) 

12 RETURN 

1991 WRITE(6, 1992) 

1992 FORMAT (’ NO. OF NEWTON ITERATION EXCEEDED LIMIT’,/, 

& ’ SOLUTION STOP’) 

STOP 

END 

SUBROUTINE OMEGAR (DG , HI , OMEGI , OMEGB , OMEGE , DC) 

C 

C CALCULATE THE BACK STRESS AT T+DT 

C 

DIMENSION OMEGB (6) , OMEGE (6) ,DC(6) , OMEGI (6) 

C 

Q1 = EXP(-DG) 

IF (ABS(DG).LE. l.E-10) DG=1.E-10 

IF (ABS(DG) .LE.l.E-4) Q2=Hl*(l.-.5.DG+DG*DG/6.-DG**3/12.) 

IF (ABS (DG) . GT . 1 . E-4) Q2 = Hl*(l.-Ql)/DG 
C 

DO 303 J=1 ,6 

OMEGE (J) = OMEGI ( J) +Q1 * (OMEGB ( J) -OMEGI ( J) ) +Q2*DC ( J) 

303 CONTINUE 
RETURN 
END 

SUBROUTINE KAPPAR (D J , H2 , Z2 , AKO , AKB , AKE , DR , MOD) 
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c 

C CALCULATE THE DRAG STRESS AT T+DT 

C 

IF (MOD. EQ. 3) GO TO 10 

AKE = AKO 

RETURN 

10 Q1 = EXP(-DJ) 

IF (ABS(DJ) .LE.l.E-10) DJ=1.E-10 

IF (ABS(DJ) .LE.l.E-4) Q2=H2*Z2* (1- .5*DJ+DJ*DJ/6. -DJ^3/12.) 
IF (ABS (D J) . GT . 1 . E-4) Q2 = H2*Z2* (1 .-Q1)/DJ 
C 

AKE = AKO+ (AKB-AKO) *Q1+Q2*DR 

RETURN 

END 

SUBROUTINE SIGMAR (DQ , KELTYP , DET , AMU , ALAM , DC , OMEGB , OMEGE , 

1 SIGB,SIGE,DVOL) 

C 

C CALCULATE STRESS AT T+DT 

C 

DIMENSION DET (6) ,DC(6) , OMEGB (6) , OMEGE (6) , SIGB (6) , SIGE (6) 

C 

PRESB = (SIGB(l) +SIGB(2) +SIGB(3)) /3 . 

Q3 = EXP(-DQ) 

IF (ABS (DQ) . LE . 1 . E-3) Q4 = 1 . -^♦DQ+DQ+DQ/6. -DQ^3/12. 

IF (ABS (DQ) . GT , 1 . E-3) Q4 = (l.-Q3)/DQ 

IF (KELTYP. EQ. 2) DET(3)=(2.*AMU^DC(3)-ALAM*(DET(1) + DET(2)))/ 
1 (ALAM+2. *AMU) 

DVOL = DET(l) + DET (2) + DET (3) 

PRESE = PRESB + (ALAM+2 . +AMU/3 . ) ♦DVOL 
C 

DO 702 J=1 , 6 
ALPHA = 1. 

IF(J, GT. 3) ALPHAS. 

DOM = OMEGE (J) - OMEGB(J) 

SIGE(J) = 2 . ★OMEGE ( J) /3 . +ALPHA*PRESE +Q3* (SIGB(J) - 
*2 . ♦OMEGB ( J) /3 . -ALPHA+PRESB) +Q4^ (2 . ♦AMU+DET ( J) - 
♦ALPHA *2 . *AMU*DV0L/3 . -2 . ♦DOM/S . ) 

702 CONTINUE 
RETURN 
END 

SUBROUTINE DELR (NIT , DVOL , ALAM , AMU , DET , SIGB , SIGE , 

1 OMEGE, AKE, AN, DT, DR) 

C 

C COMPUTE DR FOR NIT >= 2 

C 

DIMENSION DC (6) ,DET(6) , SIGB (6) , SIGE (6) , OMEGE (6) , WORK (6) 
^♦♦♦♦SECOND INVARIANT FUNCTION 

SINV(A,B,C,D,E,F)=(A*A+B*B+C*C+2.*(D*D+E*E+F*F))*2./3. 

<?♦♦ ♦♦♦HYPERBOLIC FUNCTION 
C 

DO 200 J=1 , 6 
ALPHA = 1.0 
IF(J.GT.3) ALPHA = 0. 

DC(J) = (ALPHA^ ALAM^DVOL+2 . ♦AMU^DET ( J) -SIGE ( J) +SIGB ( J) ) / 

♦(2. ♦AMU) 

200 CONTINUE 
C 

DR = SINV (DC (1), DC (2), DC (3), DC (4), DC (5), DC (6)) 

DR = SQRT(DR) 

C 
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IF(DR.LE.l.E-lO) DR=1 .E-10 

RETURN 

END 

SUBROUTINE INVER2(F,BUP) 

THIS SUBROUTINE SOLVES TWO SIMULTANEOUS EQUATIONS 

DIMENSION F (2 , 3) , FIN (2,2), BUP (3) 

FIN(1 , 1) = F(2,2) 

FIN(1 ,2) = -F (1,2) 

FIN(2, 1) = -F (2,1) 

FIN (2, 2) = F(1 , 1) 


COMPUTE THE DETERMINANT OF F 

FDET = F(1,1)*F(2,2)-F(1,2)*F(2,1) 

DO 100 1=1,2 

BUP (I) = ( -FIN (1 , 1) *F (1 , 3) -FIN (1 , 2) *F (2,3)) /FDET 
100 CONTINUE 
BUP (3) = 0.0 
RETURN 
END 

SUBROUTINE INVER3(F,BUP) 

THIS SUBROUTINE SOLVES THREE SIMULTANEOUS EQUATIONS 

DIMENSION F (3 , 4) , FIN (3 , 3) , BUP (3) 

FIN(1,1) = F(2,2)*F(3,3) - F(2,3)*F(3,2) 

FIN (1,2) = -(F (2,1)*F(3,3) - F(2,3)*F(3,1)) 

FIN(1 ,3) = F(2,1)*F(3,2) - F(2,2)*F(3,1) 

FIN(2, 1) = -(F(1,2)*F(3,3) - F(1,3)*F(3,2)) 

FIN (2, 2) = F(1,1)*F(3,3) - F(1,3)*F(3,1) 

FIN(2,3) = — (F(1,1)*F(3,2) - F(1,2)*F(3,1)) 

FIN(3, 1) = F(1,2)*F(2,3) - F(1,3)*F(2,2) 

FIN (3, 2) = -(F(1,1)*F(2,3) - F(1,3)*F(2,1)) 

FIN (3, 3) = F(1,1)*F(2,2) - F(1,2)*F(2,1) 

COMPUTE THE DETERMINANT OF F 

FDET = F ( 1 , 1 ) *FIN (1 , 1 ) +F (1 , 2) *FIN (1 , 2) +F (1 , 3) *FIN (1 , 3) 

DO 100 1=1,3 

BUP (I) = (-FIN ( 1 , I) *F ( 1 , 4) -FIN (2,I)*F(2,4)-FIN(3,I)*F(3,4)) 
k /FDET 

100 CONTINUE 
RETURN 
END 

SUBROUTINE EVALF 1 (RDOT , K , AN , AM , AN3 , AN6 , AMU , DT , F , OMEGE , 

1 AKE, AKO,DQ,DG) 

EVALUATE FI, F2 FOR WALKER’S MODEL 

DIMENSION F (2 , 3) , OMEGE (6) 
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C*****SECOND INVARIANT FUNCTION 

SINV (A , B , C , D , E , F) = (A*A+B*B+C*C+2 . * (D*D+E*E+F*F) ) *2 . /3 . 

POW = l.-l./AN 

F(1,K) = DQ-(3.*AMU*DT/AJCE)*RD0T**P0W 

ART=SINV (OMEGE (1) , 0MEGE(2) ,0MEGE(3) ,OMEGE(4) ,0MEGE(5) ,OMEGE(6)) 
IF(ART.LE.l.E-lO) ART=1.E-10 
POW = 0,5* (AM-1 . ) 

F(2,K) = DG- (AN3*RD0T+AN6*ART**P0W) *DT 
RETURN 
END 
C 

SUBROUTINE EVALF2 (RDOT , K , AN , A2 , A3 , AMU , DT , F , OMEGE , 

1 AKE,AKO,DQ,DG) 

C 

C EVALUATE FI, F2 FOR KSR’S MODEL 

C 

DIMENSION F (2, 3) , OMEGE (6) 

C*****SECOND INVARIANT FUNCTION 

SINV (A , B , C , D , E , F) = (A*A+B*B+C*C+2 . * (D*D+E*E+F*F) ) *2 . /3 . 

POW = l.-l./AN 

F(1,K) = DQ- (3. *AMU*DT/AKE) *RDOT**POW 

ART=SINV (OMEGE (1 ), OMEGE (2) , OMEGE (3) .OMEGE (4) .OMEGE (5) , OMEGE (6)) 
IF(ART.LE.l.E-lO) ART=1.E-10 
Q1 = EXP(A3*ART) -1 . 

F(2,K) = DG-A2*SQRT(ART)*Q1*DT 
RETURN 
END 
C 

SUBROUTINE EVALF3 (RDOT , K , AN , BP , HI , H2 , A 1 , A2 , Z2 , AMU , DT , F , OMEGE , 

1 AKE , AKO , DQ , DG , D J) 

C 

C EVALUATE FI, F2 A F3 FOR MILLER’S MODEL 

C 

DIMENSION F (3 , 4) , OMEGE (6) 

C*****SECOND INVARIANT FUNCTION 

SINV (A , B , C , D , E , F) = (A*A+B*B+C*C+2 . * (D*D+E*E+F*F) ) *2 . /3 . 
C*+***HYPERBOLIC SINE FUNCTION 

SH(X) = 0.5*(EXP(X)-1./EXP(X)) 

C*****HYPERBOLIC INVERSE SINE FUNCTION 
SHIV(Y) = ALOG (Y+SQRT (Y*Y+1) ) 

C 

POW = 1 . /AN 

TEMP = (RDOT /BP) **POW 

F(1,K) = DQ-3 . *AMU*DT/AKE*RDOT/SHIV (TEMP) **0.6667 
C 

ART=SINV (OMEGE (1) ,0MEGE(2) ,0MEGE(3) ,0MEGE(4) , OMEGE (5) , OMEGE (6) ) 
ART = SQRT(ART) 

AART = ART*A1 

IF(ART.LE.l.E-lO) ART=1.E-10 
IF (AART. LE. 1 .E-5) DGNEW=H1*BP*AART**AN/ART*DT 
IF (AART . GT . 1 . E-5) DGNEW=H1 *BP* (SH (AART) ) **AN/ART*DT 
F(2,K) = DG-DGNEW 
C 

T1 = H2*A2/A1*AKE**3 

T2 = H2*ART 

TP = A2*AKE**3 

T3 = H2*Z2*BP*(SH(TP))**AN 

ADIF = AKE-AKO 

IF (ABS(ADIF) .LE.l.E-6) GO TO 140 
DJNEW=( (T1-T2) *RD0T+T3) *DT/ADIF 
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GO TO 150 
140 DJNEW = 0. 

150 F(3,K) = DJ-DJNEW 
RETURN 
END 
C 

SUBROUTINE C0NST1 (DEG , EE , ANU , AKO , ANIN , AM , AN1 , AN2 , AN3 , AN4 , 

1 AN5 , AN6 , AN7 , OMEGZ , AN , ALAM , AMU , Cl , C2 , C3 , C4 , C5) 

C 

C THIS SUBROUTINE IS CALLED BY HYPELA TO CALCULATE ALL OF THE 
C TEMPERATURE DEPENDENT MATERIAL CONSTANTS FOR WALKER’S MODEL 
C 

DIMENSION TABT(6) ,EET(6) ,ANUT(6) ,AK1T(6) ,ANINT(6) ,AMT(6) ,AN1T(6) 
DIMENSION AN2T(6) ,AN3T(6) ,AN4T(6) ,AN5T(6) ,AN6T(6) ,AN7T(6) 
DIMENSION OMEGZT (6) 

DATA TABT/800., 1000., 1200., 1400., 1600., 1800./ 

DATA EET/26 . E6 , 24 . E6 , 24 . E6 , 22 . 6E6 , 18 . 6E6 , 13 . 2E6/ 

C DATA EET/26 . E6 , 24 . E6 , 23 . 4E6 , 21 . 8E6 , 19 . 6E6 , 16 . 8E6 / 

C DATA EET/26 . E6 , 24 . E6 , 23 . 4E6 , 22 . 5E6 , 21 . 6E6 , 20 . 7E6/ 

DATA ANUT/O. 322, 0.328, 0.334, 0.339, 0.345, 0.351/ 

DATA AK1T/50931. ,75631. ,95631. ,251886. ,91505. ,59292./ 

DATA ANINT/ . 059 , . 059 , . 079 , . 244 , . 195 , . 223/ 

DATA AMT/1 .158,1.158,1.158,1.158,1.158,1. 158/ 

DATA ANlT/0.,0.,0.,0.,0.,0./ 

DATA AN2T/30 . E7 , 6 . 0E7 , 1 . 5E7 , 2 . E7 , 5 . E6 , 1 . E6/ 

DATA AN3T/8000. ,1000. ,781.2,1178.6,672.6,312.5/ 

DATA AN4T/0. ,0. ,0. ,0. ,0. ,0./ 

DATA AN5T/0. ,0. ,0. ,0. ,0. ,0./ 

DATA AN6T/0 ,,0.,0.,0.,8. 977E-4 , 2 . 733E-3/ 

DATA AN7T/0. ,0. ,0. ,0. ,0. ,0./ 

DATA OMEGZT/O. ,0. ,-2000. ,-2000. , -1434. ,-1200./ 

NTP=6 

NTPM1=NTP-1 
TDIF=TABT(2) -TABT(l) 

L1=DEG 

L2=TABT(1)-TDIF 

L3=TDIF 

IT= (L1-L2) /L3 

IF(IT.LT. 1)IT=1 

IF (IT . GT . NTPM1) IT=NTPM1 

FAC= (DEG-TABT (IT) ) /TDIF 

EE= (EET (IT+1) -EET (IT) ) +FAC+EET (IT) 

ANU=(ANUT(IT+1)-ANUT(IT)) *FAC+ANUT(IT) 

AK0=(AK1T(IT+1) -AKIT(IT) ) *FAC+AK1T(IT) 

ANIN= (ANINT (IT+ 1 ) -ANINT (IT) ) +FAC+ANINT (IT) 

AM= (AMT (IT+ 1 ) - AMT (IT) ) +FAC+AMT (IT) 

AN1= (ANlT (IT+1) -AN1T (IT) ) +FAC+AN1T (IT) 

AN2= (AN2T (IT+1 ) - AN2T (IT) ) +FAC+AN2T (IT) 

AN3= (AN3T (IT+ 1 ) -AN3T (IT) ) +FAC+AN3T (IT) 

AN4= (AN4T (IT+ 1 ) -AN4T (IT) ) +FAC+AN4T (IT) 

AN5= (AN5T (IT+ 1 ) -AN5T (IT) ) +FAC+AN5T (IT) 

AN6= (AN6T (IT+ 1 ) -AN6T (IT) ) +FAC+AN6T (IT) 

AN7= (AN7T (IT+1) -AN7T (IT) ) +FAC+AN7T (IT) 

OMEGZ= (OMEGZT (IT+1) -OMEGZT (IT) ) *FAC+ OMEGZT (IT) 

AN=1 . /ANIN 

ALAM=EE*ANU/ ((1.-2. *ANU) * (1 . +ANU) ) 

AMU=(1 . -2 . * ANU) +ALAM/ (2 . *ANU) 

Cl=2. » AMU » ALAM/ (ALAM+2 . *AMU) 

C2=4 . +AMU* (ALAM+AMU) / (ALAM+2 . +AMU) 

C3=2. *AMU 
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C4=AMU 
C5=ALAM 
RETURN 
END 

SUBROUTINE C0NST2 (DEG , EE , ANU , AKO , ANIN , A1 , A2 , A3 , A4 , 

1 A5,AN, ALAM, AMU,C1 ,C2,C3,C4,C5) 

THIS SUBROUTINE IS CALLED BY HYPELA TO CALCULATE ALL OF THE 
TEMPERATURE DEPENDENT MATERIAL CONSTANTS FOR KSR’S MODEL 

DIMENSION TABT(6) ,EET(6) ,ANUT(6) ,AK1T(6) ,ANINT(6) 

DIMENSION A1T(6) ,A2T(6) ,A3T(6) ,A4T(6) ,A5T(6) 

DATA TABT/800. ,1000. ,1200. ,1400. ,1600. ,1800./ 

DATA EET/26 . E6 , 24 . E6 , 24 . E6 , 22 . 6E6 , 18 . 6E6 , 13 . 2E6 / 

C DATA EET/26 . E6 , 24 . E6 , 23 . 4E6 , 21 . 8E6 , 19 . 6E6 , 16 . 8E6/ 

C DATA EET/26 .E6,24 .E6,23 .4E6.22.5E6, 21 . 6E6,20.7E6/ 

DATA ANUT/O. 322, 0.328, 0.334, 0.339, 0.345, 0.351/ 

DATA AK1T/50931 . ,75631. ,95631. ,251886. ,91505. ,59292./ 

DATA ANINT/ . 059 , . 059 , . 079 , . 244 , . 195 , . 223/ 

DATA AlT/3 . E8 , 6 . E7 , 1 . 5E7 , 2 . E7 , 5 . E6 , 1 . E6/ 

DATA A2T/. 59,. 00179,. 66, 1.54, 14. 96, 243./ 

DATA A3T/1 .E-12, 1 .E-12, 1 .E-12, 1 .E-12, 1 -E-12, 1 .E-12/ 

DATA A4T/0. ,0. ,0. ,0. ,0. ,0./ 

DATA A5T/0. ,0. ,0. ,0. ,0. ,0./ 

NTP=6 

NTPM1=NTP-1 
TDIF=TABT(2) -TABT(l) 

L1=DEG 

L2=TABT(1) -TDIF 

L3=TDIF 

IT= (Ll-L2)/L3 

IF(IT.LT. 1)IT=1 

IF (IT . GT . NTPM1) IT=NTPM1 

F AC= (DEG- T ABT (IT) ) /TDIF 

EE= (EET (IT+1) -EET (IT) ) +FAC+EET (IT) 

ANU= (ANUT (IT+1) - ANUT (IT) ) *FAC+ ANUT (IT) 

AKO=(AKlT(IT+l) -AK1T (IT) ) +FAC+AK1T (IT) 

ANIN= (ANINT (IT+1) -ANINT (IT) ) *FAC+ ANINT (IT) 

Al= (AIT (IT+ 1) -AIT (IT) ) +FAC+A1T (IT) 

A2= (A2T (IT+1) -A2T (IT) ) +FAC+A2T (IT) 

A3=(A3T(IT+1) -A3T(IT) ) *FAC+A3T(IT) 

A4= (A4T (IT+1) -A4T (IT) ) *FAC+A4T (IT) 

A5= (A5T (IT+1) -A5T (IT) ) *FAC+A5T (IT) 

AN=1 ./ANIN 

ALAM=EE*ANU/ ((1.-2. *ANU) * (1 . +ANU) ) 

AMU= (1 . -2 . *ANU) *ALAM/ (2 . *ANU) 

Cl=2 . +AMU+ALAM/ (ALAM+2 . * AMU) 

C2=4 . *AMU* (ALAM+AMU) / (ALAM+2 . *AMU) 

C3=2 . *AMU 

C4=AMU 

C5=ALAM 

RETURN 

END 

SUBROUTINE CONST3(DEG,EE,ANU,AKO,AN,A1,A2,H1,H2,Z2,BP, 

1 ALAM , AMU , Cl , C2 , C3 , C4 , C5) 

C THIS SUBROUTINE IS CALLED BY HYPELA TO CALCULATE ALL OF THE 
C TEMPERATURE DEPENDENT MATERIAL CONSTANTS FOR MILLER’S MODEL 
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c 


DIMENSION TABT (6) , EET (6) , ANUT (6) 

DATA TABT/800. ,1000. ,1200. ,1400. ,1600. ,1800./ 

DATA EET /26 . E6 , 24 . E6 , 24 . E6 , 22 . 6E6 , 18 . 6E6 , 13 . 2E6/ 

DATA ANUT/O. 322, 0.328, 0.334, 0.339, 0.345, 0.351/ 

NTP=6 

NTPM1=NTP-1 
TDIF=TABT (2) -TABT (1) 

L1=DEG 

L2=TABT(1) -TDIF 

L3=TDIF 

IT=(L1-L2)/L3 

IF(IT.LT. 1)IT=1 

IF (IT . GT . NTPM1) IT=NTPM1 

FAC= (DEG-TABT (IT) ) /TDIF 

EE= (EET (IT+1) -EET (IT) ) *FAC+EET (IT) 

ANU= (ANUT (IT+1) -ANUT (IT) ) +FAC+ANUT (IT) 
ALAM=EE*ANU/ ((1.-2. *ANU) * (1 . +ANU) ) 

AMU=(1 . -2 . *ANU) +ALAM/ (2 . * ANU) 

Cl=2 . +AMU+ALAM/ (ALAM+2 . *AMU) 

C2=4 . * AMU* (ALAM+AMU) / (ALAM+2 . *AMU) 

C3=2. *AMU 

C4=AMU 

fiK=AI,AU 


AKO=8000 . 

AN= 1.598 
B = 1.0293E14 
HI = 1.0E7 
A1 = 9.305E-4 
H2 = 100. 

Z2 = 50000. 

A2 = 5.9425E-12 
QS = 104600. 

TM = 1588. 

TMP6 = 0.6*TM 

TK = (DEG-32 . ) *5 . /9 . + 273. 
FI = — QS / ( 1 .9859+TK) 

F3 = -QS/(.6*1.9859*TM) 

F2 = F3* (ALOG( . 6+TM/TK) +1 . ) 

IF (TK.LT.TMP6) THP=EXP(F2) 

IF (TK.GE.TMP6) THP=EXP(F1) 

BP = B*THP 

RETURN 

END 
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TABLE I. - COMPARISON OF CPU TIMES FOR WALKER'S, 
KSR'S, AND MILLER'S MODELS 
[Temperature = 1600 °F and strain rate - 3.87xlO -3 
sec -1 . ] 


Number 
of time 
steps 


CPU time 

sec 

SAFE scheme 

UVAE scheme 

ERROR1 

1x1 0 — 4 

ERROR1 
lxl 0 — 5 

ETOLB = IxlO" 2 
CTOL = 5xl0- 3 

Walker' s model 

80 

7.9 


7.4 

40 

5.8 

■T 

4 

20 

7.3 

'14 

2 

KSR's model 

80 

6 


• 7 

40 

4 

■HI ■ 

4 

20 

3 

■ 

2.4 

Miller's model 

80 

7 

22 

17 

40 

5 

21 

11 

20 

4 

+ 23 

4.5 


t Convergence is not satisfied. Fixed 
subincrement Is employed. 


TABLE II. - COMPARISON OF CPU TIME FOR 56 TIME STEPS 


Model 

CPU time, sec 


SAFE scheme 

UVAE scheme 


ERROR1 

ERR0R1 

ETOLB =0.01 


(0.0001) 

(0.00001) 

CTOL = 0.005 

Walker's 

3.91 

mmm 

4.12 

KSR's 

3.14 

HI 

4.05 
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FIG. 1 FLOW CHART OF GLOBAL- INCREMENTAL ITERATION PROCEDURE IN NON- 
LINEAR FINITE ELEfENT ANALYSIS, BASED ON INITIAL STRAIN METHOD. 
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Figure 2. - Flow chart of new uniformly valid asymptotic integration 
scheme at iocaJ level. 
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FIG. 3 HYSTERESIS LOOP PREDICTIONS OF WALKER'S MODEL 
FOR HASTELLOY-X; TEMPERATURE = 1600 °F; STRAIN RATE = 
3.87x!0~ 3 /SEC; AND NUMBER OF TIME STEPS = 80. 
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FIG. 5 HYSTERSIS LOOP PREDICTIONS OF WALKER'S MODEL FOR 
HASTELLOY-X; TEMPERATURE = 1600 °F; STRAIN RATE = 
3.87x!0‘ 3 /sec; AND NUMBER OF Tift STEPS = 20. 
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FIG. *1 HYSTERESIS LOOP PREDICTIONS OF WALKER'S MODEL 
FOR HASTELLOY-X; TEMPERATURE = 1600 °F; STRAIN RATE = 
3.87x10 3 /SEC; AND NUMBER OF TIME STEPS = 40. 
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FIG. 6 HYSTERESIS LOOP PREDICTIONS OF KSR MODEL FOR 
HASTELLOY-X; TEMPERATURE = 1600 °F; STRAIN RATE = 
3,87x!0 -3 /SEC; NUMBER OF Tift STEPS = 40. 
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FIG. 9 THER NOME CHAN I CAL FATIGUE LOOP PREDICTIONS FROM 
WALKER'S DIFFERENTIAL FORM WITH SAFE INTEGRATION 
SCHEME. 
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FIG. 1! THERMOMECHANICAL FATIGUE LOOP PREDICTIONS FROM 
KSR'S DIFFERENTIAL FORM WITH SAFE INTEGRATION SCHEME. 
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FIG. 10 THERMOMECHANICAL FATIGUE LOOP PREDICTION S FROM 
WALKER'S INTEGRAL FORM WITH THE PROPOSED UVAE INTE- 
GRATION SCHEME. 
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FIG. 12 THERMOMECHANICAL FATIGUE LOOP PREDICTIO NS FRO M 
KSR'S INTEGRAL FORM WITH PROPOSED UVAE INTEGRATION 
SCHEME. 
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(a) STACKED RING LOUVER LINER. 



(b) CLOSE-UP or LOUVERS A, 5, AND 6. 
FIG, 13 ANNULAR COMBUSTOR LINER. 



FIG. 14 POWER HISTORY FOR THERMAL CYCLE: COOLANT FLOW 
RATE, 5.5/SEC; COOLANT FLOW TEMPERATURE, 600 °F. 
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FIG. 15 CYCLIC SURFACE LINER TEMPERATURES AT TWO LOCATIONS 
ON LOUVER 5. 
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FIG. 16 COMBUSTOR LINER FINITE ELEMENT MODEL: THREE- FIG. 17 STRESS-STRAIN PREDICTION AT SEAM WELD (ELE- 

DIMENFIONAL SOLID ELEMENT WITH 546 ELEMENTS AND 1274 MENT 94, INTEGRATION PT. 1). 
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FIG. 18 STRESS-STRAIN PREDICTION AT KNUCKLE (ELEMENT 
416, INTEGRATION PT. 1). 
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